Possum_report

Исходные данные, цели и задачи

Нашей задачей было построение модели для вычисления возраста опопссумов по их морфологическим данным. Исходя из литературного поиска, мы предполагаем, что используемый набор данных был получен группой австралийских ученых в октябре-ноябре 1993 года [1, 2]. Для данного исследований были отловлены и измерены взрослые особи Trichosurus caninus, горного короткоухого поссума.

Статистическая обработка проводилась при помощи R (версия).

Загрузка необходимых пакетов производилась с помощью функции, обнаруживающей уже установленные пакеты и определяющей необходимость их загрузки:

package_installer <- function(package){
  if (!require(package, character.only=T, quietly=T)) {
    install.packages(package)
    library(package, character.only=T)
  }
}
lapply(c(library(ggplot2),
         library(dplyr),
         library(ggcorrplot),
         library(psych),
         library(stats),
         library(car),
         library(DAAG),
         library(cowplot),
         library(bslib),
         library(gridExtra)), package_installer)

В R данные являются составной частью пакета DAAG. Данные содержат следующие переменные (site — одна из семи локаций отлова поссумов: Cambarville, Bellbird, Whian Whian, Byrangery, Conondale, Allyn River и Bulburin соответственно; Pop – одна из трех популяций поссумов: Victoria, New South Wales и Queensland; sex — факторная переменная с уровнями f - самки, m - самцы; age — возраст; hdlngth — длина головы; skull — ширина черепа; totlngth — общая длина; taill — длина хвоста; footlgth — длина стопы; earconch — длина ушной раковины; eye — расстояние от медиального угла глазной щели до латерального угла глазной щели правого глаза; chest — обхват груди (в см); belly — обхват живота (в см).

Были произведены перекодировка факторных переменных и удаление переменной “Pop” (по требованию заказчика).

possum <- DAAG::possum[-3]
possum$sex <- factor(possum$sex, levels = c("m", "f"), labels = c("male", 
                                                                  "female"))

possum$site <- factor(possum$site, levels = 1:7, labels = c("Cambarville", 
                                                              "Bellbird", 
                                                              "Whian Whian", 
                                                              "Byrangery", 
                                                              "Conondale", 
                                                              "Allyn River", 
                                                              "Bulburin"))

Разведочный анализ данных

Обзорный анализ

Краткая сводка по данным:

str(possum)

## 'data.frame':    104 obs. of  13 variables:
##  $ case    : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ site    : Factor w/ 7 levels "Cambarville",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sex     : Factor w/ 2 levels "male","female": 1 2 2 2 2 2 1 2 2 2 ...
##  $ age     : num  8 6 6 6 2 1 2 6 9 6 ...
##  $ hdlngth : num  94.1 92.5 94 93.2 91.5 93.1 95.3 94.8 93.4 91.8 ...
##  $ skullw  : num  60.4 57.6 60 57.1 56.3 54.8 58.2 57.6 56.3 58 ...
##  $ totlngth: num  89 91.5 95.5 92 85.5 90.5 89.5 91 91.5 89.5 ...
##  $ taill   : num  36 36.5 39 38 36 35.5 36 37 37 37.5 ...
##  $ footlgth: num  74.5 72.5 75.4 76.1 71 73.2 71.5 72.7 72.4 70.9 ...
##  $ earconch: num  54.5 51.2 51.9 52.2 53.2 53.6 52 53.9 52.9 53.4 ...
##  $ eye     : num  15.2 16 15.5 15.2 15.1 14.2 14.2 14.5 15.5 14.4 ...
##  $ chest   : num  28 28.5 30 28 28.5 30 30 29 28 27.5 ...
##  $ belly   : num  36 33 34 34 33 32 34.5 34 33 32 ...

summary(possum)

##       case                 site        sex          age           hdlngth      
##  Min.   :  1.00   Cambarville:33   male  :61   Min.   :1.000   Min.   : 82.50  
##  1st Qu.: 26.75   Bellbird   :13   female:43   1st Qu.:2.250   1st Qu.: 90.67  
##  Median : 52.50   Whian Whian: 7               Median :3.000   Median : 92.80  
##  Mean   : 52.50   Byrangery  : 7               Mean   :3.833   Mean   : 92.60  
##  3rd Qu.: 78.25   Conondale  :13               3rd Qu.:5.000   3rd Qu.: 94.72  
##  Max.   :104.00   Allyn River:13               Max.   :9.000   Max.   :103.10  
##                   Bulburin   :18               NA's   :2                       
##      skullw         totlngth         taill          footlgth    
##  Min.   :50.00   Min.   :75.00   Min.   :32.00   Min.   :60.30  
##  1st Qu.:54.98   1st Qu.:84.00   1st Qu.:35.88   1st Qu.:64.60  
##  Median :56.35   Median :88.00   Median :37.00   Median :68.00  
##  Mean   :56.88   Mean   :87.09   Mean   :37.01   Mean   :68.46  
##  3rd Qu.:58.10   3rd Qu.:90.00   3rd Qu.:38.00   3rd Qu.:72.50  
##  Max.   :68.60   Max.   :96.50   Max.   :43.00   Max.   :77.90  
##                                                  NA's   :1      
##     earconch          eye            chest          belly      
##  Min.   :40.30   Min.   :12.80   Min.   :22.0   Min.   :25.00  
##  1st Qu.:44.80   1st Qu.:14.40   1st Qu.:25.5   1st Qu.:31.00  
##  Median :46.80   Median :14.90   Median :27.0   Median :32.50  
##  Mean   :48.13   Mean   :15.05   Mean   :27.0   Mean   :32.59  
##  3rd Qu.:52.00   3rd Qu.:15.72   3rd Qu.:28.0   3rd Qu.:34.12  
##  Max.   :56.20   Max.   :17.80   Max.   :32.0   Max.   :40.00  
## 

Для удобства дальнейшей работы с данными числовые переменные были выделены в отдельный сабсет.

possum_numerical <- possum[4:13]

Подсчет числа NA:

NA_count <- sum(is.na(possum))

NA - 3 наблюдения из 104 - менее 3%, они могут быть удалены:

possum <- na.omit(na.omit(possum))
possum_numerical <- na.omit(na.omit(possum_numerical))

Проверка на нормальность распределения данных

Проверка на нормальность распределения данных необходима, чтобы оценить круг методов, которые могут применять при дальнейшем анализе данных. При несоблюдении условия нормальности как по формальным критериям (тесты Шапиро-Уилка и Колмогорова-Смирнова), так и по визуальным характеристикам распределения (график плотности и квантиль-квантильный график) используются непараметрические статистические критерии. Они уступают по мощности параметрическим критериям, однако последние применимы только при нормальности распределения.

Построение графиков плотности

density <- apply(possum_numerical, 2, function(a) ggplot(possum_numerical, aes(x=a))+
  geom_density())

График плотности для интересующей переменной вызывается при помощи команды density$<variable>, где <variable> – имя переменной. Пример - density$footlgth. Дополнительную кастомизацию можно провести, работая, как с обычной переменной, хранящей график ggplot.

График плотности для возраста поссумов. Распределение скошено влево, не является нормальным.

График плотности для длины головы поссумов. Распределение близко к нормальному. Имеется небольшой субпик в левой части графика.

График плотности для ширины черепа поссумов. Распределение скошено влево, в правой части графика имеется небольшой субпик.

График плотности для длины тела поссумов. Распределение близко к нормальному, с несколько скошенной в правую сторону вершиной.

График плотности для длины хвоста поссумов. Распределение близко к нормальному, несколько бимодально.

График плотности для длины стоп поссумов. Распределение бимодально.

График плотности для длины ушной раковины поссумов. Распределение бимодально.

График плотности для ширины глазной щели поссумов. Распределение близко к нормальному.

График плотности для обхвата груди поссумов. Распределение близко к нормальному, вершина несколько скошено.

График плотности для обхвата живота поссумов. Распределение близко к нормальному.

Построение квантиль-квантильных графиков

Квантиль-квантильный график для интересующей переменной вызывается при помощи команды qqplots$<variable>, где <variable> – имя переменной. Пример - qqplots$earconch. Дополнительную кастомизацию можно провести, работая, как с обычной переменной, хранящей график ggplot.

Квантиль-квантильный график для возраста поссумов. Квантили указывают на распределение, отличающиеся от нормального.

Квантиль-квантильный график для длины головы поссумов. Квантили распределены нормально.

Квантиль-квантильный график для ширины головы поссумов. Распределение квантилей отличается от нормального.

Квантиль-квантильный график для длины тела поссумов. Распределение квантилей близко к нормальному.

Квантиль-квантильный график для длины хвоста поссумов. Распределние квантилей визуально несколько отличается от нормального.

Квантиль-квантильный график для длины стоп поссумов. Распределение квантилей отличается от нормального.

Квантиль-квантильный график для длины ушной раковины поссумов. Распределение квантилей отличается от нормального.

Квантиль-квантильный график для ширины глазничной щели поссумов. Распределение квантилей близко к нормальному.

Квантиль-квантильный график для обхвата груди поссумов. Распределение квантилей близко к нормальному.

Квантиль-квантильный график для обхвата живота поссумов. Распределение квантилей близко к нормальному.

Формальные тесты

Следующим шагом проверки на нормальность стало использование теста Шапиро-Уилка для указанных данных.

normality <- sapply(possum_numerical, function(a) shapiro.test(a))
normality

##           age                           hdlngth                      
## statistic 0.9374365                     0.9813663                    
## p.value   0.0001246421                  0.1648769                    
## method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a"                           "a"                          
##           skullw                        totlngth                     
## statistic 0.9381717                     0.9848761                    
## p.value   0.0001380455                  0.3045355                    
## method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a"                           "a"                          
##           taill                         footlgth                     
## statistic 0.9830628                     0.9518154                    
## p.value   0.2228095                     0.001023365                  
## method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a"                           "a"                          
##           earconch                      eye                          
## statistic 0.9171524                     0.9777459                    
## p.value   9.083728e-06                  0.08554447                   
## method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a"                           "a"                          
##           chest                         belly                        
## statistic 0.9850221                     0.9888213                    
## p.value   0.3121142                     0.5633678                    
## method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "a"                           "a"

По результатам анализа графиков и теста Шапиро-Уилка можно сделать вывод, что длина головы, общая длина, длина хвоста, ширина глазной щели, обхват груди и обхват живота распределены нормально (p-value теста Шапиро-Уилка больше 0.05, графики плотности и квантилей близки к нормальным), а возраст, ширина черепа, длина стоп и длина ушной раковины - ненормально.

Работа с выбросами

Боксплоты по количественным переменным:

boxplot(possum$age, possum$hdlngth, possum$skullw, possum$totlngth, 
        possum$taill, possum$footlgth, possum$earconch, possum$eye,
        possum$chest, possum$belly,
        names = c('age',
                  'hdlngth', 
                  'skullw', 
                  'totlngth',
                  'taill', 
                  'footlgth', 
                  'earconch', 
                  'eye', 
                  'chest', 
                  'belly'))

В переменных hdlngth, skullw, taill, eye, chest, belly отмечены выбросы в количестве, не допускающем их удаление. Далее можно проанализировать их и сделать предположения о биологическом смысле некоторых из них.

Корреляционный анализ

Чтобы при подборе модели избежать возможных проблем с коллинеарностью, был проведен корреляционный анализ данных. Была составлена матрица корреляций всех колчиественных переменных. Корреляции с p-value меньше 0.05 не считались статистически значимыми и помечались на графике крестом. Поправка на множественные сравнения не проводилась, так как результаты анализа влияли лишь на ужесточение требование к модели, а не на непосредственные выдвижение и проверку гипотез.

pleiades_mat <- corr.test(possum_numerical)
temperature_plot <- ggcorrplot(pleiades_mat$r, type = "lower", outline.col = "white", lab=TRUE, method = "circle", p.mat = pleiades_mat$p)
temperature_plot

Мы наблюдали сильную положительную связь между длиной и шириной черепа, а также между длиной стоп и ушных раковин. Мы обнаружили положительную связь средней силы между обхватом поясницы и обхватом груди, обхватом поясницы и общей длиной тела, обхватом поясницы и длиной головы, обхватом груди и длиной тела, обхватом груди и шириной черепа, обхватом груди и высотой черепа, общей длиной тела и высотой черепа. Мы обнаружили умеренную отрицательную связь между длиной хвоста и длиной ушной раковины. Также мы наблюдали умеренную положительную связь между обхватом поясницы и длиной стоп, обхватом пояницы и шириной черепа, обхватом поясницы и возрастом, обхватом груди и длиной стоп, обхватом груди и возрастом, шириной глазной щели и шириной черепа, шириной глазной щели и длиной головы, длиной стоп и общей длиной тела, длиной стоп и длиной головы, возрастом и длиной головы.

Большинство указанных корреляций было ожидаемо по биологическим причинам (пропорцианольность строения тела и его отдельных частей). Наше внимание привлекла корреляция между длиной ушной раковины и длиной стоп.

Попарные диаграммы рассеивания

Для дальнейшего анализа были построены попарные диаграммы рассеивания. На тех из них, которые содержали переменные footlhth и earconch обнаруживались два кластера. Цвет точек в первом наборе диаграмм обозначал пол, во втором - локацию поссума.

deviant_apply_sex <- apply(possum_numerical, 2, function(b) apply(possum_numerical, 2, function(a) ggplot(possum)+
        geom_point(aes(x=a, y=b, col=sex))))

deviant_apply_site <- apply(possum_numerical, 2, function(b) apply(possum_numerical, 2, function(a) ggplot(possum)+
        geom_point(aes(x=a, y=b, col=site))))

Диаграмма рассеивания для интересующих переменных вызывается при помощи команды density$<variable_y>$<variable_x>, где <variable_y> – имя переменной, которая будет располагаться по оси ординат, а <variable_x> – имя переменной, которая будет располагаться по оси абсцисс. Пример - density$footlgth. Дополнительную кастомизацию можно провести, работая, как с обычной переменной, хранящей график ggplot.

Для экономии места приведена будет наиболее показательная диаграмма рассеивания между переменными footlhth и earconch, с указанием локации измерения поссума.

Заметно, что поссумы с более длинными стопами и ушными раковинами принадлежат к популяции, обитающей в первом и втором локусе, а с более маленькими – в остальных пяти. Таким образом, можно предположить, что мы имеем дело с двумя подвидами поссумов, отличающихся по морфометрическим характеристикам и ареалам обитания.

Анализ распределения внутри подгрупп

Для более подного разведочного анализа было решено оценить распределения переменных при разбиении их на подгруппы по полу и локации.

possum$locus <- factor(ifelse(possum$site %in% c("Cambarville", "Bellbird"), "First_locus", "Second_locus"))
possum_numerical$locus <- possum$locus

Добавленная переменная соответствует гипотетическим подвидам, обитающим в одном из двух локусов.

Распределение морфометрических данных у поссумов различных полов

Распределения асимметричны, скошены влево.

Распределения в целом выглядят как нормальные

Заметная ассимметрия (скошены влево)

Распределение немного скошено вправо в случае самок, имеется два пика у самцов

Распределения не очень симметричны, в целом выглядят как нормальные

Распределения сильно отличаются от нормальных: имеют по два пика, форму, сильно отличную от нормального распределения

Распределения сильно отличаются от нормальных: имеют по два пика, форму, сильно отличную от нормального распределения

Распределения не очень симметричны, в целом выглядят как нормальные

У распределений немного искажена форма, в целом выглядят приближенно к нормальным

Распределения близки к нормальным. У поссумов женского пола имеется малый пик в правой части графика. Было высказано предположение о том, что это беременные самки.

Распределение морфометрических данных у поссумов предполагаемых подвидов

Распределения скошены влево.

Распределения в целом выглядят как нормальные Распределения близки к нормальным.

Распределение у поссумов из второго локуса близко к нормальному, у поссумов первого локуса несколько скошено вправо.

Распределения нормальны, заметно различие между двумя популяциями.

Распределения нормальны, заметно различие между двумя популяциями.

Распределения нормальны, заметно различие между двумя популяциями.

Распределение в первому локусе близко к нормальному, во втором имеет утолщенный хвост справа.

Распределения асимметричны.

Распределение в первом локусе близко к нормальному, во втором имеет утолщенные хвосты.

Результаты EDA

При статистическом анализе морфологических данных мы обнаружили выраженное бимодальное распределение по некоторым показателям (размер ушной раковины, стопы и длина хвоста). Для объяснения данной гетерогенности популяции был проведен литературный поиск, который выявил, что спустя 9 лет после сбора исследуемых данных о Trichosurus caninus, в 2002 году, было предложено разделить этот вид на северный (сохранил название Trichosurus caninus) и южный (T. cunninghami) подвиды, вследствие существенных морфологических отличий (при этом генетических оснований для разделения на отдельные виды оказалось недостаточно)[2]. То есть, данные о горном короткоухом поссуме, собранные до 2002 года, объединяют в себе сведения о двух морфологически отличных подвидах.

На рисунке 1 продемонстрированы морфологические различия этих двух видов [2]. Согласно данным австралийских коллег, у южных популяций: – более крупная ушная раковина, – значительно (P < 0,001) более длинная стопа и – значительно (P < 0,001) более короткий хвост, чем у особей из северных популяций. Животных можно четко различить с помощью индекса, основанного на этих трех морфологических показателях. Эти данные полностью согласуются с результатами нашего исследования.

Кроме того разброс значений по некоторым показателям может быть вызван примесью особей кистехвостых поссумов (Trichosurus vulpecula), поскольку они имеют схожую морфологию и пересечение ареала обитания с Trichosurus caninus.

Рисунок 1. Ключевые различия в морфологии T. cuninghami и T. caninus (по D. B. Lindenmayer et al.)

Нами также была обнаружена незначительная бимодальность в распределении величины обхвата живота: 4 наблюдения образуют вторую моду в области больших значений. Учитывая, что эти данные описывают женские особи (среди мужских особей второй моды у данного распределения не обнаружено), мы предположили, что они относятся к беременным / вынашивающим самкам. Отсутствие больших значений по этому показателю среди юных самок также согласуется с данными о возрасте половой зрелости, который наступает на третьем году жизни.

Анализ литературных данных о беременности горных поссумов показал, что внутриутробная беременность Trichosurus caninus длится чуть больше двух недель и приводит к рождению недоношенного плода очень малого размера, что вероятно, если и приводит, то к незначительному увеличению живота.

После рождения детёныш перебирается в сумку матери, где безвылазно проводит следующие 175–200 дней. Поскольку детёныш на протяжении этого срока является неотъёмным от матери, замеры живота кормящих самок либо а) проводятся с учётом младенца, либо б) отбрасываются (если дизайн эксперимента позволяет отобрать особи среди некормящих самок).

На следующем этапе вскармливания младенец больше времени проводит на спине у матери, но при этом продолжает посещать сумку до полного отъёма от груди в возрасте от 240 до 270 дней. То есть животы самок и на этом сроке выкармливания могут быть существенно увеличены за счет детеныша

Участвовали ли кормящие самки в данном исследовании не уточняется, однако, согласно данным другого исследования (проводившегося в тот же год в тех же регионах), доля вынашивающих самок в популяции составляла около 88% в летний период и 60% к концу осени (Таблица 1, [3]).

Учитывая такой высокий процент среди женской популяции, разумно предположить, что кормящие самки не были отброшены при проведении замеров. Однако вторая мода в распределении величины живота образована намного меньшей долей замеров, что не согласуется с ожидаемым высоким процентом вынашивающих самок, и может быть вызвана другими обстоятельствами: наличием нескольких детенышей в либо вариабельностью их размеров.

Таблица 1. Доля самок с детенышами (по Viggers K.L. с изменениями)

Апрель	Июнь	Август	Ноябрь
0.25(2/8)	0.88 (7/8)	0.88 (7/8)	0.60 (6/10)

Средняя продолжительности жизни у обоих подвидов Trichosurus caninus составляет 8-12 лет, при этом встречаются особи возрастом 17 лет [1]. Поссумы этого вида считаются долгожителями среди сородичей (к примеру, упомянутый выше кистехвостый Trichosurus vulpecula живёт в среднем 7-10 лет). Однако в представленных данных медианный возраст особи составил всего 3 года, а максимальный – 9 лет.

Поскольку сбор данных производился путём отлова животных в естественной среде обитания (а не в контролируемом непрерывном наблюдении), их возраст доподлинно неизвестен и был оценен косвенным способом – по степени стираемости зубов. В отличие от других измерений, полученных прямым способом (замеры частей тела).

При таком методе оценки может возникать дополнительная погрешность, вызванная травмами челюсти или генетическими дефектами эмали или пищеварения и др. Это может негативно сказаться на построении предсказательной модели.

Подбор модели

Простая линейная модель

Модели по одному предиктору:

mod <- lm(age~hdlngth, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ hdlngth, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6815 -1.3815 -0.2424  1.1479  5.0761 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -12.80894    4.79274  -2.673 0.008804 ** 
## hdlngth       0.17934    0.05165   3.472 0.000766 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.817 on 99 degrees of freedom
## Multiple R-squared:  0.1086, Adjusted R-squared:  0.09957 
## F-statistic: 12.06 on 1 and 99 DF,  p-value: 0.0007661

mod <- lm(age~skullw, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ skullw, data = possum)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.715 -1.441 -0.300  1.347  5.295 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -6.21855    3.39152  -1.834   0.0697 . 
## skullw       0.17627    0.05945   2.965   0.0038 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.845 on 99 degrees of freedom
## Multiple R-squared:  0.08155,    Adjusted R-squared:  0.07227 
## F-statistic:  8.79 on 1 and 99 DF,  p-value: 0.003795

mod <- lm(age~totlngth, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ totlngth, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2173 -1.4215 -0.1766  1.1705  4.9051 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -6.86308    3.86020  -1.778  0.07849 . 
## totlngth     0.12244    0.04418   2.771  0.00667 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.854 on 99 degrees of freedom
## Multiple R-squared:  0.07198,    Adjusted R-squared:  0.06261 
## F-statistic: 7.679 on 1 and 99 DF,  p-value: 0.006673

mod <- lm(age~taill, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ taill, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0495 -1.2322 -0.5825  1.4175  5.1840 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.50415    3.59572  -0.140    0.889
## taill        0.11676    0.09692   1.205    0.231
## 
## Residual standard error: 1.911 on 99 degrees of freedom
## Multiple R-squared:  0.01445,    Adjusted R-squared:  0.004494 
## F-statistic: 1.451 on 1 and 99 DF,  p-value: 0.2312

mod <- lm(age~footlgth, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ footlgth, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0847 -1.5317 -0.5043  1.3643  4.9591 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.07640    2.96524   0.026    0.979
## footlgth     0.05476    0.04326   1.266    0.209
## 
## Residual standard error: 1.909 on 99 degrees of freedom
## Multiple R-squared:  0.01592,    Adjusted R-squared:  0.005984 
## F-statistic: 1.602 on 1 and 99 DF,  p-value: 0.2086

mod <- lm(age~earconch, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ earconch, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9926 -1.6770 -0.6708  1.2980  5.0793 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.31802    2.28480   1.015    0.313
## earconch     0.03124    0.04730   0.660    0.510
## 
## Residual standard error: 1.921 on 99 degrees of freedom
## Multiple R-squared:  0.004387,   Adjusted R-squared:  -0.00567 
## F-statistic: 0.4362 on 1 and 99 DF,  p-value: 0.5105

mod <- lm(age~eye, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ eye, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2201 -1.4650 -0.3812  1.1994  5.2413 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -2.4912     2.6684  -0.934   0.3528  
## eye           0.4195     0.1769   2.372   0.0196 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.872 on 99 degrees of freedom
## Multiple R-squared:  0.05376,    Adjusted R-squared:  0.0442 
## F-statistic: 5.624 on 1 and 99 DF,  p-value: 0.01965

mod <- lm(age~chest, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ chest, data = possum)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7539 -1.2776 -0.1663  0.8811  4.8811 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.77202    2.43571  -1.959 0.052903 .  
## chest        0.31753    0.08975   3.538 0.000616 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.814 on 99 degrees of freedom
## Multiple R-squared:  0.1122, Adjusted R-squared:  0.1033 
## F-statistic: 12.52 on 1 and 99 DF,  p-value: 0.0006156

mod <- lm(age~belly, data=possum)
summary(mod)

## 
## Call:
## lm(formula = age ~ belly, data = possum)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.913 -1.407 -0.420  1.340  5.087 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.44667    2.15541  -2.063  0.04173 *  
## belly        0.25333    0.06581   3.849  0.00021 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.795 on 99 degrees of freedom
## Multiple R-squared:  0.1302, Adjusted R-squared:  0.1214 
## F-statistic: 14.82 on 1 and 99 DF,  p-value: 0.00021

В качестве единственных предикторов значимо предсказывают возраст переменные chest (***), belly (***), hdlngth (***), skullw (**), totlngth (**), eye (*).

Предсказательная сила простых линейных моделей оказалсь неудовлетворительной, поэтому была предпринята попытка построить множественную линейную модель для предсказания возраст поссумов.

Множественная линейная модель

Стандартизируем модель, чтобы получить информацию о том, какие из возможных предикторов оказывают наиболшее влияние на полную модель.

possum_scale <- as.data.frame(sapply(possum_numerical[-11], scale))
mod_scale <- lm(age ~ ., data = possum_scale)
summary(mod_scale)

## 
## Call:
## lm(formula = age ~ ., data = possum_scale)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.64661 -0.62064 -0.09882  0.49549  2.60641 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.804e-16  9.364e-02   0.000    1.000
## hdlngth      8.912e-02  1.648e-01   0.541    0.590
## skullw       3.445e-02  1.437e-01   0.240    0.811
## totlngth     4.928e-02  1.814e-01   0.272    0.787
## taill       -3.262e-04  1.479e-01  -0.002    0.998
## footlgth    -2.512e-01  1.916e-01  -1.311    0.193
## earconch     2.081e-01  1.808e-01   1.151    0.253
## eye          1.395e-01  1.047e-01   1.333    0.186
## chest        1.595e-01  1.495e-01   1.067    0.289
## belly        2.047e-01  1.273e-01   1.608    0.111
## 
## Residual standard error: 0.9411 on 91 degrees of freedom
## Multiple R-squared:  0.1941, Adjusted R-squared:  0.1144 
## F-statistic: 2.436 on 9 and 91 DF,  p-value: 0.01574

Заметно, что сильнее всего на модель влияют belly, earconch и footlgth.

Удаляем коллинеарные факторы посредством функции вздутия модели.

multiple_mod1 <- lm(age ~., data=possum_numerical[-11])
summary(multiple_mod1)

## 
## Call:
## lm(formula = age ~ ., data = possum_numerical[-11])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1535 -1.1886 -0.1893  0.9489  4.9917 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -1.369e+01  6.501e+00  -2.106    0.038 *
## hdlngth      4.851e-02  8.969e-02   0.541    0.590  
## skullw       2.127e-02  8.870e-02   0.240    0.811  
## totlngth     2.249e-02  8.280e-02   0.272    0.787  
## taill       -3.169e-04  1.436e-01  -0.002    0.998  
## footlgth    -1.090e-01  8.316e-02  -1.311    0.193  
## earconch     9.813e-02  8.528e-02   1.151    0.253  
## eye          2.524e-01  1.894e-01   1.333    0.186  
## chest        1.512e-01  1.417e-01   1.067    0.289  
## belly        1.437e-01  8.936e-02   1.608    0.111  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.802 on 91 degrees of freedom
## Multiple R-squared:  0.1941, Adjusted R-squared:  0.1144 
## F-statistic: 2.436 on 9 and 91 DF,  p-value: 0.01574

vif(multiple_mod1)

##  hdlngth   skullw totlngth    taill footlgth earconch      eye    chest 
##   3.0664   2.3316   3.7174   2.4696   4.1471   3.6911   1.2376   2.5238 
##    belly 
##   1.8292

multiple_mod2 <- update(multiple_mod1, .~. - footlgth)
vif(multiple_mod2)

##  hdlngth   skullw totlngth    taill earconch      eye    chest    belly 
##   3.0546   2.3233   3.4749   2.4636   1.7309   1.2375   2.4756   1.8234

multiple_mod3 <- update(multiple_mod2, .~. - totlngth)
vif(multiple_mod3)

##  hdlngth   skullw    taill earconch      eye    chest    belly 
##   2.6009   2.3231   1.4069   1.3919   1.2333   2.3826   1.8229

multiple_mod4 <- update(multiple_mod3, .~. - hdlngth)
vif(multiple_mod4)

##   skullw    taill earconch      eye    chest    belly 
##   1.8038   1.3720   1.3469   1.1760   2.2964   1.7594

multiple_mod5 <- update(multiple_mod4, .~. - chest)
vif(multiple_mod5)

##   skullw    taill earconch      eye    belly 
##   1.3433   1.3716   1.2602   1.1610   1.3685

summary(multiple_mod5)

## 
## Call:
## lm(formula = age ~ skullw + taill + earconch + eye + belly, data = possum_numerical[-11])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3800 -1.2481 -0.3102  1.0757  4.8933 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -12.24491    5.72649  -2.138   0.0351 *
## skullw        0.07380    0.06695   1.102   0.2731  
## taill         0.01277    0.10646   0.120   0.9048  
## earconch      0.03260    0.04956   0.658   0.5123  
## eye           0.24818    0.18243   1.360   0.1769  
## belly         0.18645    0.07687   2.426   0.0172 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.792 on 95 degrees of freedom
## Multiple R-squared:  0.1679, Adjusted R-squared:  0.1241 
## F-statistic: 3.835 on 5 and 95 DF,  p-value: 0.003299

После этого шага имеем модель: age = -11.64 + 0.076 * skullw + 0.003 * taill + 0.024 * earconch + 0.259 * eye + 0.183 * belly

Убираем невлияющие предикторы, которые не оказывают значимого влияния на модель.

drop1(multiple_mod5, test = "F")

## Single term deletions
## 
## Model:
## age ~ skullw + taill + earconch + eye + belly
##          Df Sum of Sq    RSS    AIC F value  Pr(>F)  
## <none>                305.20 123.69                  
## skullw    1    3.9031 309.10 122.97  1.2149 0.27314  
## taill     1    0.0462 305.24 121.70  0.0144 0.90478  
## earconch  1    1.3901 306.59 122.15  0.4327 0.51226  
## eye       1    5.9459 311.14 123.64  1.8508 0.17691  
## belly     1   18.9012 324.10 127.76  5.8834 0.01717 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

multiple_mod6 <- update(multiple_mod5, .~. - taill)
drop1(multiple_mod6, test = "F")

## Single term deletions
## 
## Model:
## age ~ skullw + earconch + eye + belly
##          Df Sum of Sq    RSS    AIC F value  Pr(>F)  
## <none>                305.24 121.70                  
## skullw    1    4.0765 309.32 121.05  1.2821 0.26034  
## earconch  1    1.4417 306.69 120.18  0.4534 0.50233  
## eye       1    5.9882 311.23 121.67  1.8833 0.17316  
## belly     1   20.8271 326.07 126.37  6.5502 0.01205 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

multiple_mod7 <- update(multiple_mod6, .~. - earconch)
drop1(multiple_mod7, test = "F")

## Single term deletions
## 
## Model:
## age ~ skullw + eye + belly
##        Df Sum of Sq    RSS    AIC F value   Pr(>F)   
## <none>              306.69 120.18                    
## skullw  1    4.2471 310.93 119.57  1.3433 0.249299   
## eye     1    5.1815 311.87 119.87  1.6388 0.203540   
## belly   1   21.9771 328.66 125.17  6.9510 0.009756 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

multiple_mod8 <- update(multiple_mod7, .~. - skullw)
drop1(multiple_mod8, test = "F")

## Single term deletions
## 
## Model:
## age ~ eye + belly
##        Df Sum of Sq    RSS    AIC F value   Pr(>F)   
## <none>              310.93 119.57                    
## eye     1     8.107 319.04 120.17  2.5551 0.113159   
## belly   1    36.141 347.07 128.68 11.3909 0.001059 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(multiple_mod8)

## 
## Call:
## lm(formula = age ~ eye + belly, data = possum_numerical[-11])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9596 -1.2969 -0.4105  1.1155  4.9715 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -7.76642    2.98114  -2.605  0.01061 * 
## eye          0.27726    0.17345   1.598  0.11316   
## belly        0.22720    0.06732   3.375  0.00106 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.781 on 98 degrees of freedom
## Multiple R-squared:  0.1523, Adjusted R-squared:  0.135 
## F-statistic: 8.803 on 2 and 98 DF,  p-value: 0.0003049

На данном шаге имеем модель age = 7.76642 + 0.27726 * eye + 0.22720 * belly

Предсказательная сила модели неудовлетворительная. Дальнейшие шаги проделаны исключительно в учебных целях.

Оценка качества модели

Строим график остатков модели

multiple_mod8_diag <- data.frame(fortify(multiple_mod8), possum_numerical[, c(2:7, 9)])

gg_resid <- ggplot(data = multiple_mod8_diag, aes(x = .fitted, y = .stdresid)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  geom_smooth(method = "lm") +
  geom_hline(yintercept = 2, color = "red") +
  geom_hline(yintercept = -2, color = "red")
gg_resid

## `geom_smooth()` using formula 'y ~ x'

С моделью не все в порядке. Имеется два наблюдения за пределами 2 стандартных отклонений (влиятельные наблюдения). Ряд измерений в центральной части производят впечатление гетероскедастичности.

Строим график расстояний Кука

cook_model <- ggplot(multiple_mod8_diag, aes(x = 1:nrow(multiple_mod8_diag), y = .cooksd)) + 
  geom_bar(stat = "identity") + 
  geom_hline(yintercept = 2, color = "red")
cook_model

График расстояний Кука выглядит приемлемо.

Проверяем нормальность распределения

quantile_model <- qqPlot(multiple_mod8_diag$.stdresid)

quantile_model

## [1]  9 11

Квантильный график выглядит не критично.

res_1 <- gg_resid + aes(x = chest)
res_2 <- gg_resid + aes(x = earconch)
res_3 <- gg_resid + aes(x = taill)
res_4 <- gg_resid + aes(x = totlngth)
res_5 <- gg_resid + aes(x = footlgth)
res_6 <- gg_resid + aes(x = hdlngth)
res_7 <- gg_resid + aes(x = skullw)

grid.arrange(res_1, res_2, res_3, res_4, res_5, res_6, res_7, nrow = 2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

Несмотря на коллинеарность, возвращаем в модель chest, поскольку наблюдаем явную неучтенную зависимость.

summary(multiple_mod8)

## 
## Call:
## lm(formula = age ~ eye + belly, data = possum_numerical[-11])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9596 -1.2969 -0.4105  1.1155  4.9715 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -7.76642    2.98114  -2.605  0.01061 * 
## eye          0.27726    0.17345   1.598  0.11316   
## belly        0.22720    0.06732   3.375  0.00106 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.781 on 98 degrees of freedom
## Multiple R-squared:  0.1523, Adjusted R-squared:  0.135 
## F-statistic: 8.803 on 2 and 98 DF,  p-value: 0.0003049

multiple_mod9 <- update(multiple_mod8, .~. + chest)
summary(multiple_mod9)

## 
## Call:
## lm(formula = age ~ eye + belly + chest, data = possum_numerical[-11])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3033 -1.2045 -0.2014  0.9790  4.8330 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -10.00075    3.27037  -3.058  0.00288 **
## eye           0.28205    0.17211   1.639  0.10450   
## belly         0.14694    0.08351   1.760  0.08164 . 
## chest         0.17668    0.11036   1.601  0.11263   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.767 on 97 degrees of freedom
## Multiple R-squared:  0.1741, Adjusted R-squared:  0.1486 
## F-statistic: 6.817 on 3 and 97 DF,  p-value: 0.000323

Несмотря на снижение значимость отдельных предикторов, общий R-квадрат несколько возрастает.

Строим график остатков для новой модели

multiple_mod9_diag <- data.frame(fortify(multiple_mod9), possum_numerical[, c(2:6, 9)])

gg_resid2 <- ggplot(data = multiple_mod9_diag, aes(x = .fitted, y = .stdresid)) + 
  geom_point() + 
  geom_hline(yintercept = 0) +
  geom_smooth(method = "lm") +
  geom_hline(yintercept = 2, color = "red") +
  geom_hline(yintercept = -2, color = "red")
gg_resid2

## `geom_smooth()` using formula 'y ~ x'

Распределение остатков стало еще хуже. За границей двух стандартных отклонений стало больше влиятельных наблюдений, а паттерн стал еще более воронкообразным.

Строим график Кука для новой модели

cook_model <- ggplot(multiple_mod9_diag, aes(x = 1:nrow(multiple_mod9_diag), y = .cooksd)) + 
  geom_bar(stat = "identity") + 
  geom_hline(yintercept = 2, color = "red")
cook_model

График расстояний Кука выглядит приемлемо.

Строим квантильный график для новой модели.

quantile_model <- qqPlot(multiple_mod9_diag$.stdresid)

quantile_model

## [1]  9 11

Квантильный график выглядит некритично.

Низкая предсказательная сила (скорректированный R-квадрат на уровне 0.15) и неудовлетворительное распределение остатков делают модель, по нашему мнению, непригодной для использования. Нижеследующий раздел (предсказания модели) выполянется в учебных целях, чтобы продемонстрировать, какой подход мы бы применяли, если бы модель была пригодна для дальнейшей работы.

Предсказания модели

MyData <- data.frame(
  belly = seq(min(possum_numerical$belly), max(possum_numerical$belly), length.out = 100),
  chest = mean(possum_numerical$chest),
  eye = mean(possum_numerical$eye))

Predictions <- predict(multiple_mod9, newdata = MyData,  interval = 'confidence')
MyData <- data.frame(MyData, Predictions)

Pl_predict <- ggplot(MyData, aes(x = belly, y = fit)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
  geom_line() + 
  ggtitle("Множественная модель")
Pl_predict

Список использованной литературы

[1] Lindenmayer, D. B., Viggers, K. L., Cunningham, R. B., and Donnelly, C. F. (1995). Morphological variation among populations of the mountain brushtail possum, Trichosurus caninus Ogilby (Phalangeridae, Marsupialia). Australian Journal of Zoology 43, 449–458.

[2] Lindenmayer, D. B., Dubach, J., and Viggers, K. L. (2002). Geographic dimorphism in the mountain brushtail possum (Trichosurus caninus): the case for a new species Australian Journal of Zoology 50, 369–393.

[3] Viggers KL, Lindenmayer DB, Cunningham RB, Donnelly CF. The effects of parasites on a wild population of the Mountain Brushtail Possum (Trichosurus caninus) in south-eastern Australia. Int J Parasitol. 1998 May;28(5):747-55. doi: 10.1016/s0020-7519(98)00022-8. PMID: 9650054.